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1. INTRODUCTION 


The numerical solution of supersonic flow problems using the full 
potential equation has become an attractive and promising alternative to 
solving either Euler equations or the linearized potential flow equations. 

The full potential equation retains most of the nonlinear features of the flow 
field, such as shocks, that the linearized potential, e.g., the more popular 
linear panel methods, inherently neglect, while having the simplicity of a 
single variable irrotational solution. Primitive variables, entropy 
singularities, and CFL stability conditions tend to complicate Euler equation 
solvers*. 

The current approach to supersonic flows was first established by 
Grossman^ for the conical flow problem using a nonconservative form of the 
potential equation. The conical flow problem reduces the full potential 
equation to an equation which in the crossflow plane (i.e., transverse plane 
normal to conical rays) contains all of the salient features of the two- 
dimensional transonic full potential equation. In Ref. 2, it was found that 
because of the type dependent, or mixed elMptic/hyperbolic nature of the 
crossflow, transonic techniques such as those developed by Jameson^ could be 
used to determine numerical solutions. The conical flow problem was extended 
to include three-dimensional flow by Grossman and Siclari^ using a fully 
implicit marching technique where each marching step requires an implicit 
crossflow solution. These schemes have all used Successive-Line-Over- 
Relaxation (SLOR) as their basis for numerical solution of the nonconservative 

T c 

full potential equation. For the nonconical problem, it was found by Siclari 3 
that accuracy could be enhanced by isentropical ly fitting the bow shock and 
numerical efficiency optimized by a judicious selection of the sweep 
direction. 

As the procedure for solving the full potential equation for supersonic 
problems matured, other investigators have established similar methods. 

Shankar et al^, using a semi-implicit marching technique with a density 
linearization of the conservative full potential equation, has also shown 
success. Comparisons of Shankar's conservative approach with the present 
nonconservative formulation in Ref. 6 show remarkably excellent agreement 
considering the conservative versus nonconservative treatments. The semi- 
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implicit formulation of Ref. 6 or 7 requires some CFL constraints, unlike the 
present fully implicit formulation which has no CFL constraints. Sritharan^ 
and Bradley et. al^ have also developed conservative formulations for solving 
the conical problem. Comparisons with the present approach that are shown in 
Reference 8 exhibit excellent agreement except for the typical 
conservative/nonconservative disparities that occur in the vicinity of the 
captured crossflow shock. 

In the present work, two basic topics are studied with the aim of 
broadening the applicability and usefulness of the present method encompassed 
within the computer code NCOREL (NONConical RELaxation, see Ref. 18) for the 
treatment of supersonic flow problems. The first topic is that of computing 
efficiency. Accelerated schemes exist and are in current use for transonic 
flow problems. One such scheme is the approximate factorization (AF) 
method. This study will develop and apply an AF scheme to the supersonic flow 
problem. The second topic that will be addressed will be the computation of 
wake flows. The proper modeling of wake flows is important for multi- 
component configurations such as wing-body and multiple lifting surfaces where 
the wake of one lifting surface will have a pronounced effect on a downstream 
body or other lifting surface. This is an interim report of a larger study. 

A final report will cover these topics and other topics presently under 
investigation. 
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2. APPROXIMATE FACTORIZATION SCHEMES 


Numerical techniques are investigated that show promise for accelerating 
convergence in comparison with the standard SLOR methods. The primary 
candidate for this is the alternating-direction-implicit (ADI)^ or, as it is 
more commonly referred to in its application to nonlinear transonic flows, 
approximate factorization (AF) schemes. These AF schemes have been applied 
successfully to a variety of transonic flow problems. Initially, AF schemes 
were applied to the Transonic Small Disturbance (TSD) equation by Ballhaus, et 
al^. Holst^ successfully applied an AF2 type scheme to the conservative 
full potential equation for transonic flows. The nonconservative full 
potential equation was treated successfully by Baker^»14 f or 20 transonic 
flows and should be applicable to the nonconservative full potential 
supersonic/transonic crossflow problem of the present study. Two basic AF 
algorithms are considered, ADI or AF1 and the AF2 scheme which splits one of 
the second derivatives into two first-order derivative operators. The latter 
scheme has reportedly been the most stable in supersonic flow regions for the 
transonic flow problem. 


2.1 BASIC FORMULATION 

The nonconservative form of the 3D full potential equation is written in 
a spherical coordinate system (w, r). The governing equation is then 

transformed via a conformal stereographic projection to (p, q, t) coordinates 
and further by a crossflow conformal Jowkowski mapping to (p, e, R) 
coordinates. In terms of a reduced potential F, where Q = v<j> + q,,, and $ = 
RF(p,e,R), the full potential equation can be written as 


< a2 - v2 > F pp + - 


■ RH [(W Z -a Z )RHh 1 . ZVIVHhjF^ ♦ f F p „ + F^ 


+ i [<W 2 -a Z )RHh 2 ♦ 2MU][h 1 F pe ♦ / F e0 ♦ 

+ ■* H (H 2 -a 2 )[l. 1 f pR */F 9R + F RR ] + ... (1) ; 
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where U, V, and W are the velocities in the 8 , p, and R directions, 
respectively, and H is the combined metric of the two mappings. In general, 
the conformal mapping in the crossflow plane leads to nonorthogonal coordinate 
derivatives if the mapping singularity is a function of r. This mapping 
dependence on r leads to mesh derivatives defined as hj = Pp and *13 = p8 r . 

The radial marching direction r, remains unchanged due to the transformations, 
or r = t = R. The details of the mappings and coefficients can be found in 
Ref. 2 and 4. 

Unlike transonic flow, the supersonic flow problem is contained within a 
finite crossflow mesh bounded by a bow shock. The bow shock may be captured 
within a prescribed outer boundary 2 * 4 or the bow shock can be fitted as the 
outer boundary 5 using the isentropic shock jump conditions. ! 

A shearing transformation is applied to Eq ( 1 ), between the body p = B(e, 
R) and outer boundary or bow shock p = C ( 8 , R), where 

X = 8 



Z = R 


which yields a final rectangular computational mesh. 

Equation (1) can be represented as the sum of a conical plus a nonconi cal 
operator, in the form 

LUij) = L c (* i#j ) + R • L NC (* 1#j ) . (3) 

The nonconi cal coefficients on the RHS of Eq ( 1 ) all have an R dependence and 
vanish identically at R = 0 for the quasi -two-dimensional conical flow 
problem. 

2.2 CONICAL FLOWS 

For conical flows, after applying the shearing transformation, Eq ( 1 ) can 
be written as 


L c Ui,j) = A i F xx + A 2 F XY + A 3 F YY + **• 


(4) 
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where 


A, = 


isizyii 

2 

p 


;2UV y ^laWly 

2 p P p 2 0 


A = (a 2 .v 2 )Y 2 - My y + li^idlY 2 
a 3 la ; R P pe ^2 0 


( 5 ) 


Equation (4) closely resembles the nonconservative form of the 2-D 
transonic flow equation. The difference is that the type dependency of the 
conical part of Eq (1) or Eq (4) is linked to the nature of the crossflow 

6 ? o 

velocity defined by Q = Ir + \r. An upwind bias in the difference equations 

c .22 
must occur when the crossflow velocity is supersonic or Q c > a . The 

crossflow velocity component V is always negative and toward the body 

surface. The U component of velocity can be positive or negative depending 

upon the geometry and angle of attack. Equation (4) has heretofore been 

2 4 

solved successfully using transonic SLOR techniques * and the rotated 

3 

difference scheme of Jameson . 

The principal part of Eq (4) can be rewritten in a rotated difference 
2 2 

format as, for Q > a 


where 


A 1 F XX + A 2 F XY + A 3 F YY = (a2 " Q c } F ss +a2p nn 


U 


ss 


v ? 2U 1 V 1 

1 C 4. 1 1 c + 1 p 

0 2 yY q 2 r XY 7? XX 

X x x 


nn 


U« 21LV, Vo 

0 2 ^ YY Q 2 r XY Q 2 r XX 

X X X 


( 6 ) 
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and 


U x = U/p , 

v 2 = V/p , 

An upwind bias is applied to the finite difference representation of the F ss 
terms, and central differences are used for the F nn terms in supersonic 
crossflow regions. 

The following sections will present an adaptation of the two basic AF 
algorithms, ADI or AF1 and AF2, to the present supersonic flow problem. 

2.2.1 AF1 Factorization 

An ADI or AF1 type factorization can be applied to the principal terms of 
Eq (4), for subsonic crossflow Q c 2 < a 2 , in the form 


UY 

V, = V Y + — i 

1 p p 

VY 

IL = UY 

2 p p 



a “ L c ( »i J 1 


(7) 


Where A i , j is the correction to the reduced potential F i j or Aj j = F i j n+1 - 
Fi,j n * w is a relaxation parameter and L c ( ^ j ) is the residual of Eq (4) at 
the nth iteration, a is an acceleration parameter which is varied in a cyclic 
fashion during the iteration process. The two first-order difference 
operators result in a second-order central difference. The first-order 
operators are defined as 


t 


*x = < 


'i+l »J " ( } i,j 


^i-l»j ’ ^ ^i ,0 


^i.j+1 


( 




1 J-l 


( 8 ) 
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The basic premise behind an approximate factorization scheme can be 
revealed if the LHS of Eq (7) is expanded and terms not resembling those of Eq 
(4) are neglected, or 



Vx . . Vy 

? + m 3 
AX^ 6 




9 i ,j 


(9) 


If a) = 1, Eq (9) would be equivalent to solving Eq (4) with the cross 

derivative evaluated using old values of the potential. Equation (9) is 

typically solved in a two-step format by defining an intermediate variable, 

Gi .• , where 
* 



G*? + ^ = cuoL ( d>? • ) 
1 ,J C M ,J ' 



( 10 ) 


Equation (10) represents two tri diagonal systems of equations involving 
differences only in the computational X or Y direction. Equation (10) must be 
modified in regions of supersonic crossflow to include the proper upwind 
bias. The following form of the supersonic crossflow factorization is 
essentially identical to the AF1 scheme of Baker^*^ equation. Hence, for 



K l A lu 



K 2 A 1u 




= cuoL( j) 


( 11 ) 


a - A 


3c 


AY 2 


- A 



.n+1 r n+l 
A. . = G . . 
1 ,J 1,J 


The central (subscript c) and upwind coefficients (subscript u) are given 
by the rotated difference scheme Eq (6) as 
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'lc 


5 ? 


(a‘ 


'lu 


Qq) 2 

— Uj 


, = U 2 
3c 7 C 


(a 2 - Q 2 ) 2 

Su = X s - V 1 


( 12 ) 


The first factor allows for the U component of velocity to be positive or 
negative, where if 

U > 0, K x = 1, K 2 = 0 

(12a) 

U < 0, = 0, K 2 = 1. 

As mentioned earlier, the V velocity is always negative for supersonic 
crossflow and, hence, only forward upwind differences occur in the second 
factor. 

In general , the first factor involves a pentadi agonal matrix and the 
second factor a quadradi agonal matrix. As suggested by Baker 13 , these 
differences can be replaced by 


Vy ■ < >?!)-« >?:)-< + < 


,n+l 




= { ) n+1 - ( ) n+1 - ( ) n + ( ) n 

i 1 i+l,j 1 i+1 ,j 1 i+2,j 


(13) 


Vy ~ ( } i J " ( l 1 J+l " ( 'ij+l + { ^i ,j+2 ‘ 


This reduces the set of Eq (11) to the following tri diagonal form, for 


2 2 
Q 2 > a 2 , 


< , . - W T A lu 

1C ^' 


t 


(K 1 S X + K 2 S X») 


( 14 ) 


t t 

Vy . . U Y \ n+1 _ r n+l 

3c aY 2 3u aY 2 / 1 » J 1 ,J 


where, for U > 0, I s = 1 and U < 0, I s = -1 
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2.2.2 AF2 Factorization 


11 12 

The AF2 algorithm has been found by several investi gators 11 to 

generally be the more stable AF factorization for transonic flows. In the AF2 

factorization, one of the second derivatives is split between the two 

factors. Following a similar AF2 factorization used by Baker , for subsonic 
2 2 

crossflow, Q c < a , the AF2 algorithm becomes 



AX 2 


a + 



“ L c (*i.j ) 


(15) 


This form of factorization is thought to be more stable since the Y 
operator in the first factor yields a <f>Yt term, unlike the ^ term of the AF1 
scheme 14 . The term H m in the two factors accounts for the transformation 
derivatives of the mapped and sheared mesh. As illustrated by Catherall 1 ^, 
the proper splitting of these transformation derivatives can yield optimum 
convergence, and neglecting these derivatives can considerably degrade 
convergence. The coefficient A 3 in Eq (4) contains the shearing transforma- 
tion derivatives Y p and Y 0 /p. The mapped olane (p,0) velocities, V and U 
respectively, contain the metric H(p,0) of the two mappings. Hence, a 
suitable form for the factor^ H m might be 



AY Y p AY 1 fv V 
AX H or AX H lY p ‘ p 1 


(16) 


The first-order forward difference operator on Y is placed in the first factor 
since this term does not switch for supersonic crossflow conditions. For 
supersonic crossflow, the AF2 factorization is modified to include an upwind 
bias and a tri diagonal form as 



^('Vx* K 2 dS X>) G iJ 


awL ( 4>? . ) 

c vr i ,r 


a + A 3c H m aT + A 3u H m aT 


A? +1 = G? +1 . 
i.J 


(17) 
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Since the forward Y operator is included in the first factor, the sweep is 
constrained to be toward the body surface or decreasing J. 

An alternate factorization to the AF2 scheme described above would be to 
split the X-derivative. Unlike transonic flow, the difficulty in splitting 
the X-derivative in the supersonic problem is that the U velocity in the 
supersonic crossflow region can be either positive or negative. This occurs 
primarily at the captured bow shock. In transonic flow, the X or U velocity 
direction is much like the Y direction of the present problem in that a 
negative supersonic U velocity is unlikely to occur. Hence, there is no first 
order X operator that, in general , does not switch. One could propose a 
scheme where the factorization is set up for U <0, and if u > 0, the upwind 
coefficient is either neglected or a shift operator is imposed. This scheme 
was found to be unstable or would not work for this problem. 

Other AF2 factorizations were considered, including the AF3 factorization 
of Baker 13 where both coefficients are brought into the first factor. This 
would seem to be a candidate for a faster scheme since the differential 
operators would not act on the coefficients and lead to spurious terms. 

Baker , in fact, has reported his AF3 scheme to be considerably faster than 
either the AF1 or AF2 schemes. Unfortunately, this scheme could not be 
applied successfully to Eq (17). 

2.2.3 Boundary Conditions 

Figure 1 illustrates the conformally mapped and sheared computational 
crossflow plane domains. Symmetry conditions are imposed at e = ± w/2 or 
I = 2 and IC for the symmetric half plane problem. Hence, periodic end 
conditions apply on Y = constant lines. On X = constant lines, J = 2 
corresponds to the body surface and J = JMAX corresponds either to the outer 
boundary (BSC) or the bow shock (BSF) . In both the bow shock capture or fit 
methods, the outer boundary has the same condition that the correction to the 
potential vanish, or 


A dummy row or Y = constant line at J = 1 is used to implement the body 

i ^ flow tangency. This condition relates the values of the 

correction at J = 1 to those at J = 3, or for conical flow yields: 


10 


( 19 ) 


In this way, central derivatives can be used for Fy and the body surface 
coordinate line can be treated like any other coordinate line. 

. The order of the factors in both AF schemes were chosen so that the first 
sweep is carried out on Y = constant lines. This was chosen over the reverse 
factorization so that periodic end conditions could be imposed on the 
intermediate variable G, or 


G, . = Gq . 

1,J 3,j 

G IC+l,j = 6 IC-l,j 


( 20 ) 


which most certainly is a reasonable assumption. If the factors in the AF 
schemes are reversed, then somewhat arbitrary boundary conditions must be 
imposed on the intermediate variable G. The end conditions (18) and (19) then 
apply to the second sweep on the X = constant lines. 

2.2.4 Temporal Damping 

It has been indicated that the basic AF1 scheme may be unstable in 
supersonic regions. To allow for this possibility, the AF1 scheme has been 
generalized to include an explicit temporal damping F st (e.g., <j> s j.) term. 
Jameson's generalized AF scheme^ includes similar terms in both factors. In 
Reference 12, it was indicated that this term may also be required in the AF2 
scheme. It was found that adding this term explicitly to the first factor was 
sufficient to maintain stability for flows with captured shocks. This 
stability or temporal damping term has the format 

V 1 U 1 

“ G 6 st = " e O^T + QjT ( K 1^X + K 2^ 

where 


K 

jng 


< e < e K 
max 


e =10 
max 


(21) 
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In general, the addition of the temporal damping slows convergence. Hence, 
the form of the factor e was chosen to maximize the damping in the vicinity of 

sonic lines or across shocks so as not to cause an overall degradation in the 

convergence rate. The constant K is chosen to be as small as possible for the 
optimum convergence. 

In general, there are no restrictions on the sweep due to the velocity 
directions in an AF1 scheme. The first sweep can be from the outer boundary 
to the body surface or the reverse. In order to properly include temporal 
damping in the AF1 scheme, the first sweep must be in the direction of the 
supersonic crossflow. This requires that the Y = constant sweep must be 
towards the body or decreasing J since V < 0. 

In the supersonic problem (unlike the transonic case), it was observed 

that the AF1 scheme required little or no temporal damping except for the most 
difficult cases whereas the AF2 scheme required considerably higher values of 
e for the strong Y-shock solutions. 

2.2.5 Acceleration Parameter 

For the AF1 scheme, the maximum and minimum values of the acceleration 
parameter a were taken as 

“max ~ A max 
“min “ A min 

The coefficients A fn ^ n and A max were taken anywhere from 0.5 to 4.0 for 
all the cases computed. The convergence rate could be affected by as much as 
a factor of two by a judicious choice of these parameters. 

In the AF2 scheme, these parameters were chosen as 
“max “ A max 

“min ” A min 


1 

a5T 


1 

’ W 


1 

a5T 


1 

AY 


(23) 
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? * 7 

AX 4 AY^ 

( 22 ) 

1 1 
AX ’ AT • 
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where varied between 0.5 and 6 and A^* between 3 and 6. Typically, 
unity for A min and three for A max was sufficient for most cases. The 
theoretical value of = 1 or 2 could never be achieved, possibly due to 
the effect of the various coordinate transformations. 

For both AF schemes, the cyclic variation of a took the form 


a 


I 


“min 
amax “max 


1-1 

IMA'X-T 

for I = 1, IMAX 


(24) 


where IMAX = 3 for the AF1 scheme and IMAX = 6 for the AF2 scheme. The 
minimum number of cycles seemed to be the best choice for the AF1 scheme. The 
convergence rate of the AF2 scheme was affected insignificantly for IMAX 
between 3 and 6. Further increase in IMAX # reduced the convergence rate. In 
all of the cases computed, w^pj = 1.50 or 1.70 and u^p 2 = 1.33. Departures 
from Eq (24) did not seem to affect significantly the convergence rates. 

2.2.6 Conical Results 

Two techniques are available for the computation of supersonic conical 
flows: the bow shock capture^ (BSC) and bow shock fit^ (BSF) methods. In 
the BSC method, an outer mesh boundary p = C( 6) is prescribed, and the bow 
shock is captured within this boundary. The BSC method is a more stringent 
test of the AF schemes in that two shocks may be present in the flow, the bow 
shock (Y-shock) and an embedded crossflow shock (X-shock) as illustrated in 
Fig. 1. The bow shock is the most critical in that its position and strength 
largely determine the internal flow. The bow shock also extends around the 
entire field encompassing more points than the embedded crossflow shock. The 
BSF method fits the bow shock as the outer boundary and, hence, eliminates the 
bow shock from the internal flow calculation; if an embedded supersonic 
crossflow region is not present, the internal flow problem becomes elliptic. 
The conical convergence rate of the BSF method is largely determined by the 
implicit shock fitting procedure where the shape of the boundary is updated 
and usually underrelaxed until the isentropic shock jump condition is 
fulfilled at each shock mesh point. The conical BSF method also requires more 
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computational time per iteration because the crossflow mesh and metrics must 
be recomputed for each iteration after the bow shock shape has been updated. 

For conical flow, the BSC method is used primarily to evaluate the AF 
schemes. Figure 2 shows the effect of reversing the order of the factors in 
the AF1 scheme for a thin elliptic cone at = 2.0, a = 0°. The AF1XY scheme 

represents the order of the factors indicated in Eq (7) and (14) where 
periodic end conditions are used for the intermediate variable G. The factors 
were then reversed (AF1YX) with the first sweep occurring on X = constant 
lines. Two different boundary conditions were used: setting j = 0 as the 
end condition in the tri diagonal on the dummy row below the body; and using 
^i»l = ^i»3 as the P en °dic end condition, as is the case for the condition on 
the correction ^ . 

All three cases were run with the same a variation. As shown in Fig. 2, 
the results are quite sensitive to the boundary condition on the intermediate 
variable, and making the intermediate variable mimic the correction seems to 
be the best choice. Even with this boundary condition, the YX factorization 
does not give identical results to the XY factorization and seemed to be 
somewhat more sensitive to the a variation. Hence, the AF1XY scheme was used 
for all the computations. 

Figure 3 shows a comparison of the convergence rate of the maximum residual 
for the AF1 and SLOR schemes for a thin subsonic leading edge elliptic cone (0 C 
= 20°, <S C = 2°) on a (48 x 38) crossflow mesh at = 2 and angles of attack a = 
0° and 10°. The SLOR scheme found to be optimum for this problem in Ref. 5 is 
one that sweeps around the body on X = constant lines. The "column" SLOR scheme 
was found to be two to five times faster than the alternate SLOR scheme which 
sweeps toward the body on Y = constant lines. It is interesting to note that 
after one or two orders of magnitude reduction in the maximum residual, the SLOR 
scheme. for the supersonic freestream problem does not exhibit the typical 
slowdown in convergence rate that occurs in transonic flow problems. A break. in 
the SLOR curve occurs after one order of magnitude but remains linear for 
further reductions. It is also interesting to note that the SLOR convergence 
rates of the a = 0° and a = 10° flows are not dramatically different, both 
taking about 350 to 400 iterations, considering that the 10° case is a 
multishocked flow. At a = 10°, a strong crossflow shock develops but is 
evidently overshadowed by the convergence of the captured bow shock. The AF1 
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AF1 Scheme (BSC) Convergence Rate 
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Fig. 3 Comparison of AF1 Convergence Rate with SLOR (BSC) 
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scheme at a = 0° converges very quickly. Essentially, these flows are 
converged when the maximum residual reaches 10“^. For a = 0°, this occurs at 
about 10 iterations or when the log (RESMAX) = -2. The AF1 scheme is an order 
of magnitude faster than the SLOR scheme for a = 0°. As the angle of attack 
increases, the AF1 scheme slows down by a factor of 3 while the SLOR remains 
about the same. Overall, the AF1 scheme is at least three times faster 
iteration wise than the SLOR scheme. The relaxation factor was 1.5 for these 
cases in both schemes, and three cycles were used in the AF1 scheme. A larger 
number of cycles did not seem to enhance convergence. Figure 4 shows the 
surface pressure distributions for the elliptic cone computed in Fig. 3 at a = 
0°, 5°, and 10° and = 2.0 on a finer mesh. Both a = 5° and a = 10° have 
crossflow shocks on the leeward surface. Figure 5 shows the computed 
crossflow streamline pattern at a = 10°. As mentioned earlier, the V 
component of velocity is negative. The crossflow streamlines emanate from the 
bow shock and travel toward the body surface coalescing at the leeward and 
windward vortical singualrities. One streamline stagnates on the body and 
wets the body surface. Also shown is the extent of the embedded supersonic 
crossflow region, which terminates at the crossflow shock. 

Before the AF2 scheme was implemented, the effect of splitting the 
transformation derivatives between the two factors was studied. Figure 6 
indicates the effect of using different forms for the term H m in Eq (16). If 
the transformation derivatives are neglected in the factorization, the case 
could only be run when the a variation was increased significantly. 

Increasing the minimum value of a generally degrades the convergence rate. 

The best convergence was achieved when both the metric and the shearing 
transformation derivatives were included in the term H m - The two curves in 
Fig. 6 with H m other than unity were obtained with an alpha variation that 
diverged when H m = 1. Hence, the AF2 scheme seems to be sensitive to the 
coordinate transformations, and the convergence rate can be affected 
significantly. The form of H m is not considered to be optimum, and further 
analytical and numerical studies should be conducted to study its effect on 
the convergence rate. A nonoptimum H m may also affect the minimum values of a 
that can be used. 

Figure 7 shows a comparison of the convergence rates of the AF1 and AF2 
schemes for the elliptic cone of Fig. 3 at a = 10°. A comparable convergence 
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rate was eventually achieved with the AF2 scheme. Further investigation of 
the AF2 scheme led to more problems (see Ref. 17) and more sensitivity (i.e., 
not user friendly) to the choice of acceleration parameters. Hence, the AF2 
scheme does not seem to have any advantages over the AF1 scheme and was not 
further considered for implementation in NCOREL. 

2.2.7 Mesh Refinement 

In Ref. 14, it was observed that for transonic flow, mesh refinement 
enhanced convergence for SLOR but did not significantly accelerate convergence 
with the AF schemes. This is apparently not the situation with supersonic 
flows where the bow shock formation is a critical factor. This is 
demonstrated in Fig. 8 for the subsonic elliptic cone at a= 5°. The mesh 
refinement in the AF1 scheme yields a factor of three enhancement in 
convergence rate over the same computation on a single fine mesh. Mesh 
refinement also enhances the convergence rate of the SLOR scheme. The mesh 
refined AF1 indicates a factor of six to seven over the mesh refined SLOR 
result, which is similar to the reduction gained for a single mesh solution. 

2.2.8 Bow Shock Fit 

The AF1 scheme was also adapted to the conical bow shock fit method. In 
the BSF method, the capture of the bow shock is eliminated from the internal 
flow field. For any given outer boundary within or coincident with the true 
bow shock, the internal flow field will quickly converge. Hence, the 
convergence rate is largely governed by the iteration scheme for determining 
the bow shock position that satisfies the isentropic shock jump conditions. 
Typically, the updated bow shock boundary must be underrelaxed at each 
iteration so as not to disturb the convergence of the internal flow field 
computation. The bow shock boundary also has to be underrelaxed so as not to 
overshoot the correct bow shock position, which would cause the internal flow 
attempt to capture a portion or all of the bow shock, thus leading to 
divergence. The major advantage in applying the AF schemes to the BSF method 
is that the correct bow shock information will propagate faster to and around 
the boundary and allow greater values of the shock relaxation parameter (uj s ), 
leading to an enhanced overall convergence rate. Figures 9 and 10 illustrate 
the AF1 convergence rate of the subsonic elliptic cone in comparison to SLOR 
at a = 5° and 10°, respectively. In Fig. 9 and 10, the maximum residual 
convergence is shown. In the SLOR curves, the shock relaxation parameter was 
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0.75 and 0.50, respectively, for a = 5° and a = 10°. Higher values cause 
divergence of the flow field. It was observed that with the AF1 scheme the 
shock relaxation parameter did not require under-relaxation for these cases 
and could be taken to be equal to unity. Hence, in the AF1 scheme both the 
bow shock and internal flow field converge faster than the SLOR. As a result, 
a similar enhancement in convergence rate is obtained for the BSF method as 
compared to the earlier results indicated for the BSC method. 

2.3 NONCONICAL OR THREE-DIMENSIONAL FLOWS 

The AF1 and AF2 schemes both worked well for quasi -two-dimensional 
conical flow yielding convergence rates two to ten times faster than SLOR. 

The AF2 scheme was somewhat more sensitive because of the split Y derivative 
and the need to include the coordinate transformation derivatives in the 
factorization. Hence, for the present study, the AF1 scheme will be 
considered for the nonconi cal or three-dimensional flow problem. When the 
shearing transformation is applied to Eq (1), the principal terms of the 3-D 
full potential equation can be rewritten as 

L(<i>i,j) = (A x + B 1 )F XX -+ (A 2 + B 2 )F xy 

+ (A 3 + B 3 ) F yy + BgF zz + B g F xz (25) 


+ B 10 F yz + B n F z + 


where the B^ coefficients represent the additional nonconi cal R or 3-D 
terms. The Bj coefficients are rather complicated and are defined in Ref. 4. 

The geometry is assumed to be conical at the apex or R = 0 of the 
configuration. Marching solutions are then obtained on spherical R or Z = 
constant surfaces. The terms F^ and F^z always have upwind differences, 
whereas the F^ and F Y £ terms are smoothly switched when mixed subsonic/ 
supersonic crossflow occurs. A first-order accurate Fjz difference requires 
information at two previous planes. Initially, the AF1 scheme was applied to 
the crossflow plane XY terms of Eq (25) with the Z derivatives treated as 
forcing terms evaluated with old values of the potential. This scheme turned 
out to be slower and resulted in divergence in many cases when compared to the 
optimum SLOR, which includes all the principal Z terms in the tridiagonal 
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matrix. Hence, a factorization was sought that would maintain the XY 
crossflow convergence rate of the AF scheme and still retain the SLOR 
efficiency for the Z terms. The following AF1Z factorization was proposed for 
subsonic crossflow, < a 2 : 


(A x + B x ) 


Vx ~ B 8 ' B 9 6 Xc - B ll 


AX 




G ?J " 


° " ^ A 3 + B 3^ 
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_ r n+l 
" G i,J 


(26) 


where 6 Xc and 6 yc represent second-order, central, first-derivative operators, 
or 


6 Xc ( ) i+l,j ' ( ) i-i,j 


6 Yc = ( ) i,j+l - ( ) i,j-l * 


(27) 


The Z terms (e.g., F 22 ) do not factor because of the hyperbolic nature of the 
problem and hence are added explicity to the two factors. 


The F^^ difference contains the unknown value of the potential at the 
current station and the two known values of the potential at the two previous 
stations. Except for the nonconical coefficients of the F xx and F YY terms, 
the Z terms are added explicitly to each factor. Including the Z terms in the 
factorization was also necessary to maintain diagonal dominance and conver- 
gence. In many cases, neglecting the Z terms in the AF scheme not only slowed 
convergence but actually resulted in divergence. For supersonic crossflow, 

the factorization was modified to include the upwind Z terms or, for O 2 > a 2 

9 x c * 



<Ai u +B,> 


AX 


r~ (k iV k 2 S x’ 


B 8 B 9 B ll\ r n+l 

- aKat - srj G i,j 


a 


- A. 



(A 3U +B 3> 


$Y ^ Y 

' B 10 a7a I 



_ r n+l 
■ G 1,0 


ouoL( .) 

1 >J 
(28) 


25 


It was found that the off diagonal terms of the F X z derivative could be 
included in the subsonic crossflow region but not in the supersonic region 
leading to the following for the 6 XZ operator for U > 0, or 


*XZ ' ( *1 J,k ‘ ( l 1-lJ,k " ( ,j ,k-l + ( Vl.j.k-l 


(29) 


K-, = ( )?!} 4 . - < )? 4 U - ( > 


xz 


i+l,j,k ‘ ' 'i ,j ,k ~ 1 'i+l,j,k-l +( ) i,j,k 


where the subscript K refers to the present R station and K-l refers to the 
known values of the potential at the previous station. Hence, care must be 
taken to preserve the proper balance of new and old values of the potential. 


Because of the split nature of the governing equation and R dependence, 
the AF1Z scheme was found to converge most reliably and optimally if the 
acceleration parameter a was scaled with R, or 


“min ~ “min TflC * aT {1 + R) 
“max ~ %ax » ^7 ( 1 + R) 


(30) 


Hence, the acceleration parameter variation reduces to the conical value at R 
= 0 and increases linearly with R as the marching solutions are obtained. 
Since the nonconi cal coefficients are scaled with R, the acceleration 
parameters must also be scaled. This scaling will cause the retention of 
approximately the same convergence rate for a body of length unity versus an 
arbitrary dimensional length. 


2.3.1 Three-dimensional Applications of the AF1Z Scheme 


Although it is impossible to test the AF1Z scheme for all possible 
situations, test case computations were carried out on a variety of arbitrary 
wings and bodies. The following test cases were computed with the bow shock 
fit option. 


Figure 11 illustrates the convergence history of the first case of a 
highly swept delta wing at = 1.80, a = 0°. The wing geometry consists of a 
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Fig. 11 Squire Wing Computation, AF1Z vs SLOR (BSF) 
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parabolic centerline thickness distribution with elliptic spanwise cross 
sections. As is the case for most three-dimensional wings, the geometry 
commences near the apex with a thick (e.g., 3:1 major to minor axis ratio) 
cross section which becomes thin and eventually approaches a flat plate cross 
section at the trailing edge. Figure 11 shows the number of iterations 
required per marching step to reduce the maximum residual to 10-3 (i.e., a 
minimum of four orders of magnitude) for the SLOR and AF1Z schemes. Step 0 
refers to the conical or R = 0 solution required to start the computation. 
Fitting the bow shock as the outer boundary reduces the internal flow field 
computation to an elliptic problem if an embedded supersonic crossflow region 
does not form. For this case, the entire flow field was elliptic. Mesh 
refinement is used only at the conical station, and the marching proceeds on 
the fine mesh using the previous station solution as an initial guess. The 
iterations at step 0 reflect only the fine grid convergence. The AF1Z scheme 
shows large gains for the initial steps. At the initial steps, the geometry 
changes most rapidly. The AF1Z gains taper off as the cross section becomes 
quite thin. The SLOR scheme has some difficulty computing the latter stations 
in comparison with the AF1Z scheme. Step 30 corresponds to the centerline 
trailing edge. The calculation is taken beyond the trailing edge in order to 
compute the entire wing. The wake is assumed to be a flat plate, which in 
this case is an exact assumption. The total iterations for the run are also 

shown in Fig. 11. The SLOR computation took 2010 fine grid iterations, and 

the AF1Z scheme required only 672 iterations, producing an overall factor of 
three reduction in iterations. The actual computation time, which includes 
geometry and mesh generation, was reduced by a factor of two. 

The next case, shown in Fig. 12, is for a 67° sweep angle arrow wing at 

M» = 1.75, a = 5°. The geometry consists of a symmetrical, NACA 4% thick, 
four-digit airfoil imposed chordwise on the wing. The wake is approximated as 
a flat plate. In this calculation, the flow field very quickly becomes 
supercritical at R = 3 or STEP 4. The same trends apply for this case except 
that the SLOR scheme does not have any rise in iterations near the trailing 
edge of the wing. The interesting aspect of the AF1Z scheme is that the 
number of nonconi cal iterations required per step is relatively constant and 
independent of the geometry variation. Almost a factor of three reduction in 
iterations is again achieved by the AF1Z scheme corresponding to a factor of 
two in running time. Hence, the appearance of supercritical crossflow and a 
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crossflow shock does not seem to deteriorate the AF1Z scheme significantly. 

In addition, aft of the trailing edge, the crossflow shock merges with the 
trailing edge shock. 

Another case, shown in Fig. 13, is for a realistic, supersonic maneuver, 
demonstration wing designed with the aid of NCOREL^® and tested by NASA 
Langley. Details of this wing can be found in Ref. 19. This wing has a 
variable sweep leading and trailing edge. The leading edge planform angle 
varies from 25° to 33°. The wing also has significant twist and camber. 

Figure 13 shows the convergence histories for the two schemes at Mo = 1.62 and 
a = 14°. The wake is approximated by a flat plate extension of the wing 
spanwise camberline. For this wing, supercritical flow also appears at R = 3 
or STEP 4. Similar trends apply for this wing except that the AF1Z scheme 
gains are reduced aft of the trailing edge of the wing. The overall reduction 
in iterations is similar, being almost a factor of three and corresponding to 
a factor of two reduction in running time. 

All of the previous cases were run on a relatively fine grid (58 x 58). 
Figure 14 shows the iteration ratio of the SLOR to AF1Z scheme for several 
grids, corresponding to the DEMO wing computation of Fig. 13. The AF1Z scheme 
performs almost as well on the cruder meshes as on the finest mesh, retaining 
between a factor of two to three reduction in number of iterations. Figure 14 
does indicate that the performance gain will increase on the finer grids. 

The AF1Z scheme performs well for wing computations. The next set of 
cases will test the scheme for body computations. Figure 15 shows the 
convergence history for an axi symmetric circular arc cylinder body at M ra = 
1.60, a = 10° on a 58 x 58 mesh. For this case, the crossflow is subsonic 
along the entire length of the body and, hence, shockless. The AF1Z scheme 
performs nicely by reducing the number of iterations by a factor of four. 

Figures 16 to 18 show another set of computations for a more difficult 
body shape: an elliptic cross section (3:1 axis ratio) with a Haacke-Adams 
area distribution. Figure 16 shows the computations at = 1.60, a - 5° on a 
58 x 58 mesh. Better than a factor of four reduction in number of iterations 
is achieved for this condition. Figure 17 shows a higher Mach number case (M„ 
= 2, a = 5°) for the same body. For this condition, just under a factor of 
three reduction was achieved. Both schemes seem to have difficulty in the 
supercritical crossflow region that commences at around 20 in the marching 
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Fig. 12 Arrow Wing Computation, AF1Z vs SLOR (BSF) 
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Fig. 17 Elliptic Missile Body Computation at M^ = 2.0, a = 5°, AF1Z vs SLOR 


steps. Figure 18 shows an even more extreme and difficult condition to 
compute, M^ = 2, a = 10°. For this condition, the crossflow is supersonic at 
the conical start, with the crossflow shock becoming increasingly stronger 
toward the aft end of the body. The AF1Z scheme performs well initially, but 
both schemes begin to have difficulty towards the aft end of the body. Not 
quite a factor of two reduction in iterations is achieved. 

In summary, the AF1Z scheme seems to be reasonably reliable and can be 
made to run faster than the SLOR scheme over a wide range of configurations. 
Even for difficult bow shock fit computations, the AF1Z will run faster 
although the gains achieved are not as great as at the lower Mach number 
conditions. It must be mentioned that the effect of the AF parameters was not 
intensively explored for the previous computations and, hence, a more thorough 
study of these parameters might yield greater reductions. 

2.4 PARAMETER SELECTION 

Although it is impossible to test the AF1Z scheme for all possible 
situations, test case computations were carried out on a variety of arbitrary 
wings and bodies. The a variation that was found to work best in a variety of 
cases was 


a 

max 


-h (l + R) 

ay^ 


_ l 

“min "AX’ 


(1 + R) 


(31) 


where the cycles between and are two or three for crude to medium 

meshes (16 x 16 + 48 x 48) and four for finer meshes (58 x 58 and above). The 
exception to this rule is for more difficult bow shock fit cases at high Mach 
number or incidence. For these cases, the alpha range of Eq (31) should be 
raised, for example, in Eq (31), A„, 1n ■ U \, ax - 2; these values might be 
raised to A min = 2, A max = 4 for the crude mesh to get the bow shock fit 
computation started. In addition, the number of cycles ITMAX might be raised 
to 6 on the crude mesh. After the crude mesh is successfully computed, the 
default values can probably be used on the finer and nonconical meshes. 

The AF1Z scheme can be run with a relaxation factor (w) for the residual 
up to 1.7. Since the AF1Z scheme converges faster than the SLOR scheme, the 
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bow shock relaxation (u> s ) factors can also be increased both at the conical 
start and at the nonconi cal stations. Typically, io $ = 3.0 was used at the 
nonconical stations. The temporal artificial viscosity or damping terms 
controlled by the parameter EST are also included in the AF1Z scheme. For 
most cases, with the exception of difficult cases with strong embedded shocks, 
the AF1Z sheme can be run with the EST parameter set to zero. For cases that 
require damping, the value of EST should be set at very small values (e.g., 
-0.001 > EST > -0.10). For most cases -0.01 or -0.001 is sufficient. Nonzero 
values of EST in the AF1Z scheme degrade the convergence considerably and 
should be avoided. 
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3. WAKE FLOWS 


In general, wake flows are characterized by a discontinuity or slip in 
velocity with a continuity in pressure; In potential flows, a jump in 
potential is prescribed across the wake which accounts for the circulation or 
lift about the wing or airfoil. In two-dimensional potential flows the 
velocity is continuous, but the potential has a constant jump across the wake 
streamline. Matching pressures at the trailing edge in the wake cut also 
impose a Kutta condition and cause the wake streamline to leave the trailing 
edge in a direction along the trailing edge bisector. For the Euler equations 
the situation can be different because of entropy losses across shock waves. 

If shock waves of different strengths exist on the upper and lower surfaces of 
an airfoil a slip line with a jump in velocity will arise in the wake due to 
the differences in total pressure and entropy along the upper and lower 
surface streamlines. 

On the other hand, in three-dimensional flows, contact discontinuities 
arise in the wake for both Euler and potential flow formulations. In 
addition, in three dimensions, depending on the trailing edge geometry and 
flow conditions (e.g., cusped or finite angle), the trailing edge streamline 
will leave tangentially to one of the surfaces or at the mean angle of the 
trailing edge^. In potential flows, the slip in velocity is due to a 
discontinuity in the direction of the total velocity vector at the wake 
surface. The magnitude of the velocity vector must match on the wake to 
impose the Kutta condition and continuity of pressure along the wake 
surface. The discontinuous direction of the wake surface total velocity 
vector causes at least two of the three components of velocity to be 
discontinuous or to have a slip. For the most part, in three-dimensional 
transonic flows, the wake is treated in a similar fashion to its two- 
dimensional counterpart. The jump in potential for a particular airfoil 
section is assumed constant to downstream infinity. This approximation 
matches the longitudinal or axial velocity. A continuous velocity is also 
imposed through the wake surface. The spanwise velocity on the wake surface 
is not matched and is just proportional to the variation in circulation or 
lift and, hence, to the variation in the spanwise potential jump for each 
airfoil section. Hence, the magnitude of the total velocity vector and 
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resultant pressure is only approximately matched on the wake surface. This is 
a good approximation for wings that are not highly swept or tapered. A 
similar wake treatment is used by Shanker, et al 2 * for three-dimensional 
supersonic flows. 

Up to the present time, no attempt had been made to model the wake in 
NCOREL, and a flat plate treatment had been utilized. The successful 
application of NCOREL to three-dimensional wings and bodies suggested that a 
wake treatment might further extend its applicability to wing-body and more 
complex configurations. 

In the present formulation, the wake is also modeled approximately as a 
planar cut in three dimensions. . Here, the Kutta condition at the trailing 
edge and the continuity of pressure on the wake surface are imposed exactly. 
Therefore, as distinct from other treatments, the behavior for the jump in 
potential is not prescribed a priori, but rather computed at every point on 
the wake surface. Initially, a conical wake problem will be discussed along 
with some special aspects concerning this problem. These basic techniques are 
then readily extendable to the fully three-dimensional flow about wings and 
wing-body configurations. 

3.1 CONICAL WAKE PROBLEM 

A conical wake problem was devised to test the basic numerics required to 
solve the fully three-dimensional supersonic wake flow. In this problem, 
sketched in Fig. 19a, the wing cross section and wake surface share a common 
apex. This problem has already been considered by J.H.B. Smith^ 22 ^, where 
some general behavior was postulated theoretically and will be discussed in a 
later section. 

The wake is modeled as a planar surface cut of infinitesimal thickness. 
The wing cross section and wake cut (see Fig. 19b) are then mapped to a near 
circle (p, e, R) using conformal mappings and further sheared to a 
computational domain (X, Y, Z). 

Flow tangency is obtained by using a dummy row of grid points around the 
body whose values are obtained by imposing the vanishing of the normal 
velocity through the Fy derivative. For positive angle of attack, the flow on 
the lower surface of the wake cut is computed with the boundary condition that 
the velocity component through the wake cut is continuous. 
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The velocity components for conical flow are defined in the mapped 
coordinate system (p, e, R) as: 


u = iff F e + u » = iff (F x x e +F y Y e ) + u » 
v = ff F P + v ~ = ff F y Y P + v ~ 


( 32 ) 


W = F + w 
00 


where $ = RF(p,e) and U is circumferential and W is the radial velocity 
component. The airfoil and wake cut are defined as Y = 0 in the computational 
domain. If a variable jump in potential exists along the wake segment, in 
general, the U velocity will reflect this discontinuity or slip, and the V 
velocity (which does not have F x in its definition) will be constrained to be 
continuous. 


For V to be continuous, the condition that 



(33) 


must be imposed on the wake cut where the change in sign reflects the opposite 
sense of the mapped cylindrical coordinate system. This condition (33) 
further reduces to 

F.,Y - F Y + 

(-Tr) = - (~Tr) (34) 


si nee the correspond!' ng upper and 1 ower wake surface poi nts map to the same 
point in physical space yielding the same freestream velocity. The wake 
condition (34) is imposed by computing the value of Fy in a one-sided fashion 
on the upper wake cut and is then used as the boundary condition on the lower 
wake cut. The computational method proceeds by computing on lines X = 
constant commencing at the lower symmetry plane and sweeping around to the 
upper symmetry plane. On the lower wake cut, the governing full potential 
equation is satisfied subject to condition (33). On the wing surface, flow 
tangency is satisfied. Both conditions use a Neumann boundary condition given 
by the dummy row of potentials. On the upper wake cut, the full potential 
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equation is not satisfied. Instead, a jump in potential is assumed to exist 
at every point on the upper wake cut as 

f ;- f ; + 4 f( v- 

The upper wake cut X = constant lines are then solved as a Dirichlet problem 
once a jump in potential has been imposed. 

In summary, the X = constant lines emanating from the lower wake cut are 
solved by using a Neumann type boundary condition obtained for condition (33) 
where V w + is derived from a one-sided difference in the upper wake plane. 
Hence, communication across the wake exists without explicitly differencing 
across the wake cut, thus eliminating the necessity for interpolation. An 
equation for the jump in potential is now needed that matches the pressures 
all along the wake cut. For potential flows, it is a sufficient condition to 
match the total speed along the wake cut for the pressures to match. Equality 
of total speed can be expressed by the equation 


0(u w + - u w-> + «v-v + 8( v-v = 0 

where 


(35) 


0 = U + U 
w+ w- 


V = V + V 
w+ w- 


8 ■ M w + W w- • 

Some interesting aspects of the conical wake problem can be uncovered by 
inspection of Eq (35). At the symmetry planes = 0 and Vy^ = -V w _; Eq 
(35) then reduces to 


U 


w+ 



= ± tt/2 


which implies, given the definition (32), that 
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Hence, in the conical wake problem, for the pressures to be matched across the 
wake cut at the symmetry plane, the jump in potential must vanish. Thus, at 
the trailing edge, the wake initially reflects the jump in potential coming 
off the trailing edge of the airfoil, and eventually the jump must vanish at 
the symmetry plane in order for the pressures to match. 

In practice, Eq (35) can be implemented recursively by solving for tf in 
terms of the U, V and average velocities on the wake cut. The new value of 
the jump in potential is computed based on old values of the velocity and is 
under relaxed until the full potential equation, flow tangency, and conditions 
(33) and (35) are all satisfied. 

Figure 20a shows the conical pressure solution on an (82 x 58) mesh about 
the airfoil and wake cut at M*, = 2.0, a = 5° and 0^ = 20°, &■[•£ = 10°. A 
crossflow shock occurs on the upper surface of the airfoil. The pressures 
match smoothly at the subsonic trailing edge. A grid point is not situated 
exactly at the trailing edge. A slight expansion occurs near the symmetry 
plane along the wake cut. Figure 20b shows the corresponding potential 
distribution about the airfoil along with the jump in potential computed for 
the wake cut. The jump in potential along the wake cut is relatively linear 
except near the symmetry plane. Figure 20c shows the overall computed isobar 
pattern. Figure 20d shows the streamline pattern in the vicinity of the 
airfoil. The trailing edge streamline leaves the airfoil smoothly, and it 
appears to leave tangentially to the lower surface. The streamlines crossing 
the wake cut are discontinuous due to the slip in velocity on the wake 
surface. An interesting feature of this computation is that the wake 
circumferential V velocities are negative along a portion of the upper wake 
cut near the symmetry plane. The wake streamline terminates at the symmetry 
plane as a nodal or vortical singularity. 

Figure 21 shows another computation for a thinner airfoil section at = 
1.7, a = 5° with 0 LE = 30° and 0y E = 20°. This computation behaves in a 
similar fashion exhibiting a subsonic trailing edge behavior. For the 
trailing edge to be supersonic in the conical wake problem, the crossflow Mach 
number at the trailing edge would have to be supersonic. The author has not 
been able to find a conical situation in which a supersonic crossflow exists 
at the trailing edge. 
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Fig. 20 Conical Wake Solution for - 2.0, ot - 5 (#• c -* 20 , 9 y£ 10 ) (Sheet 1 of 2) 







It was mentioned previously that the conical wake problem has been 
discussed by J.H.B. Smith 22 . Based on conical flow theory, Smith postulated a 
crossflow streamline pattern and nodal behavior such as is sketched in Fig. 

22. Smith postulated that the wake streamline would roll up into a nodal 
point. Hence, the vorticity shed by the airfoil would be carried to 
downstream infinity via the radial velocity at the crossflow nodal or 
stagnation point. This behavior can not be obtained within the present 
formulation, but attempts at fitting the conical wake problem indicate that 
this behavior may occur. The wake could be fit at low angles of attack (e.g., 
1° or 2°) but could not be fitted at higher angles of attack, possibly due to 
the beginning of a rollup behavior. 

3.2 GENERAL THREE-DIMENSIONAL WAKES 

Unlike the conical wake problem, the jump in potential at each point on 
the wake surface emanating from a three-dimensional wing has as its origin an 
upstream trailing edge point. Hence, for a lifting wing, in the crossflow 
plane, a variable jump in potential exists along the wake surface proportional 
to the circulation about an upstream streamline section. The conical 
condition that the jump in potential must vanish at the symmetry plane no 
longer applies to the three-dimensional wing. This condition is alleviated in 
fully three-dimensional flows because the radial velocity component now 
contains R derivatives of the potential that relate the jump in potential at 
the symmetry plane to its upstream origin. The boundary conditions and 
relationships matching the pressure are the same in the three-dimensional 
problem as in the conical problem with the exception that the velocities have 
three dimensional or R derivatives included in their definition. 

In some respects, the three-dimensional problem is somewhat simpler to 
compute since a good initial guess for the jump in potential in the crossflow 
plane can be generated from the jump in potential across the wake and airfoil 
section of the previous plane. The variation of the jump in potential along 
the wake surface is not nearly as large as that which occurs in the conical 
wake problem. Some underrelaxation must be used in updating the potential 
jump but the overall rate of convergence is not affected significantly except 
when a strong shock interaction occurs at the trailing edge. 

In general, both subsonic and supersonic trailing edges can occur in the 
3-D supersonic wake flow problem. Subsonic trailing edges have local Mach 
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numbers in a plane normal to the trailing edge that are less than unity. For 
this situation, the wake streamline leaves the trailing edge smoothly. When 
the local normal trailing edge Mach number is supersonic, the airfoil surface 
streamline will either expand or pass through a shock in satisfying the 
pressure condition. Depending upon the value of the local normal Mach number 
at the upper and lower surface trailing edge points, several flow situations 
can exist. For supersonic trailing edges, these are 

o Shock - shock 

o Shock - expansion. 

A shock or expansion can also occur on one surface in combination with smooth 
subsonic behavior on the other surface. Figure 23 sketches some of the basic 
flow situations and their character in the spherical crossflow plane. The 
trailing edge shock is really a three-dimensional surface which takes a shape 
similar to that shown in Fig. 23a in the crossflow plane at zero incidence. A 
crossflow shock can also coexist on the surface of the cross section that will 
interact with the trailing edge shock(s) and/or expansion. One such 
complicated interaction is sketched in Fig. 23d from observations of computed 
crossflow plane isobar patterns. The crossflow shock and trailing edge shock 
intersect, forming two resultant shocks. The crossflow shock after 
intersection with the trailing edge shock exists in the flow on the wake cut 
( i . e - , passes through the trailing edge shock) and is attenuated below the' 
wake cut by the lower surface expansion. 

, l 

3.2.1 Symmetric Cross Sectional Geometries 

Symmetric geometries were treated initially. The meshes generated for 
these geometries by NCOREL yield corresponding grid points in the upper and 
lower surfaces of the wake cut. Hence, no interpolation is required at 
corresponding upper and lower wake points for the potential or the speed. 
Figure 24 shows a selected sample of the crossflow plane surface pressure 
solutions for a symmetric arrow wing with 70° leading edge and 45° trailing 
edge sweep angles at M,,,, = 1.75, a = 5° commencing with the spherical surface 
that cuts through the centerline trailing edge. Some of the computed wake 
pressure distributions are compared with the pressure distribu-tions obtained 
by assuming a flat impermeable plate for the wake. A crossflow shock exists 
on the upper surface of the airfoil . The upper surface pressures indicate a 
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Fig. 23 Sketch of Trailing Edge Shock Interactions (3-D) in the Crossflow Plane 
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trailing edge shock. The lower surface pressures indicate a smooth flow 
behavior or a very slight expansion. The flat plate solution indicates shocks 
on both surfaces of the airfoil. Further down the wing, the trailing edge 
shock becomes somewhat stronger, and the crossflow shock approaches the 
trailing edge of the wing. Finally, both shocks merge just as the crossflow 
shock intersects the trailing edge and form a single strong shock at the 
trailing edge. At the last station, the crossflow shock has passed through 
the trailing edge shock and exists on the wake cut. The trailing edge 
solution indicates a shock-expansion behavior at this point. The trailing 
edge shock, crossflow shock, and expansion interaction are similar to that 
sketched in Fig. 23d. 

Figure 25 indicates a similar behavior at M,,, = 1.75, a = 10° for the same 
arrow wing. The two shocks merge causing a very strong shock at the trailing 
edge, and eventually the crossflow shock passes through the trailing edge 
shock and sits on the wake cut with a shock expansion at the trailing edge. 
Naturally, this interaction has the built-in approximation of modeling the 
wake cut as a planar surface. What the effect of fitting the wake exactly 
would have on these complicated interactions is not known at this point. 

The important aspect of modeling the wake properly is to be able to 
predict the effect of the wake of a lifting surface on a fuselage or other 
lifting surfaces. To test the wake model, a series of computations were 
carried out on a set of arrow wings built and tested by NASA Langley (see Ref. 
23). Figure 26 shows a sample of the isobar patterns computed, along with the 
planform shape and centerbody for two of the four models tested at Mach 
numbers of 2.36, 2.96, and 4.63. Model 1 has a leading edge sweep of 63.4°, 
and Model 2 was highly swept with a leading edge sweep of 76°. The models 
were instrumented for pressure, and detailed pressure measurements are 
available on both the wing and centerbody at flow spanwise sections. One of 
the spanwise stations was downstream of the centerline trail in<} edge and 
serves as a test for the accuracy of the present wake model . One of i:hese 
wings was computed previously using NCOREL without a wake model, with the 
result that the body pressures could not be predicted (see Ref. 24). The 
planform isobar patterns clearly show the crossflow shock on the leeward 
surface. The! crossflow shock intersects the trailing edge wherein the 
complicated interaction of Fig. 23d takes place. Figure 27 shows a comparison 


56 





p J 
ro 


D) R = 30 


1264-035 (2/2) (T) 


C) R = 28 

Fig. 25 Sample Arrow Wing Computation at M m = 1.75, a = 10° (Sheet 2 of 2) 
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Fig. 26 Isobar Planform Pattern for NASA Arrow Wings 
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of measured and computed pressures for Model 1 at two different conditions, M^ 
= 2.36, a = 6° and M w = 2.96, a = 10°. Excellent agreement is achieved at the 
lower Mach number on both the wing and centerbody. At the higher Mach number 
and incidence, slightly lower pressures are computed in the vicinity of the 
upper surface trailing edge. Excellent agreement is achieved for the lower 
wing surface. The body pressures are in good agreement except near the upper 
surface symmetry plane. Unfortunately, the resolution of the body is quite 
poor on the mesh that is currently used because the relative size of the body 
with respect to the wing is quite small for Model 1. Comparisons for the 
highly swept Model 2 at a = 3° and 6° are shown in Fig. 28 and 29 for M«, = 

2.36 and 2.96, respectively. Better resolution of the body is obtained for 
this model because of its relatively large size in comparison to the wing. 

Good' to excellent agreement is achieved for body pressures. Good agreement is 
achieved for the lower surface of the wing. Poor correlation is achieved on 
the upper surface of the wing at the higher incidence. The higher pressure 
supercritical plateau shown by the measured data at a= 6° is usually 
indicative of leading edge flow separation and vortex formation and, hence, 
correlation with computed pressures would not be expected. 

3.2.2 Arbitrary Cross Sectional Geometries 

. In general, the wing cross-sectional geometry is not symmetric and can 
have, camber or twist associated with it. If the leading edge is dropped, the 
placement of the singularity of the conformal wing mapping will generate a 
grid that will cause a translation of the grid points in the wake cut. Hence, 
a lower wake cut grid point will not have a corresponding upper wake point. 
This internal grid generation complicates the implementation of the wake 
computation. Interpolation of the potential and speed from the lower wake 
mesh locations into the upper wake must be carried out in order to match the 
pressures and transverse velocity in the wake cut. 

Figure 30 shows a selected sample of the crossflow plane solutions for a 
cambered arrow wing at M,,, = 1.75, a = 5°. This is the same arrow wing as Fig. 
24 but with leading edge chordwise droop implemented with a parabolic chord- 
wise camber. Figure 30a shows an early station in the wake. The translation 
of the upper and lower wake cut points is evident. This becomes increasingly 
apparent farther downstream in the wake, where twice as many points appear on 
the wake cut. The interpolation scheme seems to work well for general 
geometries and has also been implemented for wing-body geometries. 
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Fig. 28 Model 2 Comparison with Experiment at Spanwise Wake Station 
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Fig. 29 Model 2 Comparison with Experiment at Spanwise Wake Station 
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Fig. 30 Cambered Arrow Wing (£2| p = 70°, J2 TP = 45°) M m = 1.75, a = 5' 
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4. CONCLUSIONS 


An approximation factorization technique has been successfully applied to 
the three-dimensional supersonic flow problem. An average overall factor of 
two reduction in computer time can be expected from the AF1Z scheme. 

Supersonic wake flows have also been computed in both conical and three- 
dimensional full potential flows by modeling the wake as a planar cut with a 
potential discontinuity. Continuity of pressure was satisfied exactly on the 
wake cut by computing the necessary potential jump at each wake cut grid 
point. The conical wake problem is unique in that the jump in potential 
vanishes in the crossflow symmetry plane. Generally, good correlation was 
achieved for fully three dimensional wake flows with the planar wake 
approximation. 
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APPENDIX A: 


NCOREL 

PHASE II VERSION 
USER'S MANUAL UPDATE 
(NASA CR-3694) 


NOTE: The following variables have been added to the NAMELIST input. Section 

VIII summarizes all of the additional namelist inputs. 


I. Approximate Factorization Scheme (AF1Z) 

An AF1 scheme has been added to NCOREL for both bow shock capture and fit 
options. This scheme can reduce the computer time by a factor of two to 
three. 


Namelist Variable 

IAF 

IAF = 0 (SLOR) 

IAF = 1 ( AF1Z) 

Conical Control Parameters 

NCYC(K) , K-l, KREF 
AFMIN(K) , " 

AFMAX(K), " 

NCYC = 3 or 4 
AFMIN = .5 
AFMAX =2.0 


Nonconi cal Control Parameters 

NCYCR 

AFMINR 

AFMAXR 


Description 

Optional switch for relaxation scheme 


Description 

NCYC is the number of cycles from the 
maximum to minimum acceleration 
parameter. 

High Mach (i.e., > M^ = 2.5) starts 
for bow shock fitting on crude grids 
may require 6-8 cycles and AFMIN =3, 
AFMAX = 6. Sensitivity mostly to 
AFMIN. 


These parameters treated similar to 
conical parameters. 


Important: The AF1Z is sensitive to the temporal damping parameter (EST). 

: Small values of EST (i.e., -.01 or -.001) should be used or 

convergence will deteriorate dramatically. 


II. Crossflow Y-Stretching 

NPOW 

NPOW = 0 EVEN 

NPOW = 1 Clustered 

toward body 


Two types of crossflow radial meshes 
are available. 

NPOW = 0, evenly spaced radial mesh; 
NPOW = 1, hyperbolic stretching which 
clusters grid points near the body 
surface. 
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III. Bow Shock Fitting Parameters 
I ENTRY 


IENTRY = 0 
IENTRY = 1 

ITSHKC 

ITSHKR 


ITSHKC = 0 to 15 
ITSHKR = 0 to 5 


The code will fit the bow shock 
isentropically or with the Rankine- 
Hugoniot shock conditions. With the 
Rankine-Hugoniot shock, the flow 
field pressures are corrected by the 
computed bow shock entropy. 

Isentropic bow shock 
Rankine-Hugoniot bow shock 

The number of initial iterations for 
which the bow shock conditions are 
held fixed before updating. 

ITSHKC, conical 
ITSHKR, nonconical 


IBOW 


I BOW = 1 


An additional option has been 
included to implicitly fit the bow 
shock in 3D. This option yields a 
fully second order bow shock 
computation. 

First order bow shock. Marching 
meshes held fixed. 


Second order implicit bow shock. 3D 
mesh position updated in a similar 
fashion to conical. 


Important: IBOW - 2 will significantly increase the computing time by a 

factor of two to three. In addition, it should only be used at 
high Mach (i.e., M„, > 2.5). Can eliminate first order marching 
instability at high Mach number using large step sizes. The 3D 
bow shock relaxation parameter (RELNC) must be reduced to near 
conical values for IBOW = 2. 

*SM00 Bow Shock Angle Smoothing. This 

option should only be used as last 
resort if bow shock instability 
arises at high Mach number. 


ISMOO = 0 
ISMOO = 1 


Can be used with IB0W=1 or IB0W=2. 

No smoothing. 

Smoothing. 
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IV. Geometry Interpolation 


ILIN 

ILIN = 0 


ILIN = 1 


Optional mapped space body shape 
interpolation. 

Linear Interpolation. Should be used 
when cross sections have abrupt 
changes. For example, unfaired wing- 
body cross sections or nearly 
rectangular body shapes. The cubic 1 
spline fitting can produce 
unpredictable results. The maximum 
number of geometry points (e.g., 99) 
should be used with ILIN = 0. 

Cubic spline fitting. 


V. Supersonic Crossflow Region 

ISUP 

ISUP =1 First order accurate treatment of 

supersonic crossflow points. 

■ >..i Slip -.2 Second order accurate treatment of 

supersonic crossflow points. 


VI. Wake Computation 

IWAKE 
IWAKE = 0 

IWAKE = 1 

I BODY 
IBODY = 0 
IBODY = 1 


Wake cut treated as impermeable flat 
plate 

Turns on wake computation. Wake 
should be modeled as smooth extension 
of trailing edge. 

Only used if IWAKE f 0 
Wing alone 

Wing-body configuration 


Note: If IBODY = 1, a subroutine BODYW must be included that yields the 

maximum width of the body at a given axial location. 


RELWK 

i 

.25 < RELWK < 1.0 


Relaxation factor for wake pressure 
match 
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VII. Other Namelist Parameters 
DM I NR 

l.E-4 < DMINR < l.E-2 
0ME6R 

1.5 < OMEGR < 1.7 


Maximum residual tolerance for 
convergence at nonconi cal stations. 

Relaxation parameter for potential 
at nonconi cal stations. 


Namelist Parameters 

IAF 


VIII. Summary of Additional 

1. Approximate Factorization 

2. Grid Control 

3. Bow Shock Fitting 

4. Geometry Interpolation 

5. Supersonic Crossflow 

6. Wake Computation 

i 

7. Other 


NCYC(K) ,AFMIN( K) ,AFMAX(K) , K=l,3 
NCYCR, AFMINR, AFMAXR 

NPOW 

IBOW 
I ENTRY 
ITSHKC 
ITSHKR 
ISMOO 

IL IN 

ISUP 

I WAKE 
IBODY 
R£LWK 

DMINR 

OMEGR 
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